Image charge dynamics in time-dependent quantum transport 
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In this work we investigate the effects of the electron-electron interaction between a molecular 
junction and the metallic leads in time-dependent quantum transport. We employ the recently 
developed embedded Kadanoff-Baym method [Phys. Rev. B 80, 115107 (2009)] and show that 
the molecule-lead interaction changes substantially the transient and steady-state transport prop- 
erties. We first show that the mean-field Hartree-Fock (HF) approximation does not capture the 
polarization effects responsible for the renormalization of the molecular levels neither in nor out of 
equilibrium. Furthermore, due to the time-local nature of the HF self-energy there exists a region in 
parameter space for which the system does not relax after the switch-on of a bias voltage. These and 
other artifacts of the HF approximation disappear when including correlations at the second-Born 
or GW levels. Both these approximations contain polarization diagrams which correctly account for 
the screening of the charged molecule. We find that by changing the molecule-lead interaction the 
ratio between the screening and relaxation time changes, an effect which must be properly taken 
into account in any realistic time-dependent simulation. Another important finding is that while in 
equilibrium the molecule-lead interaction is responsible for a reduction of the HOMO-LUMO gap 
and for a substantial redistribution of the spectral weight between the main spectral peaks and the 
induced satellite spectrum, in the biased system it can have the opposite effect, i.e., it sharpens the 
spectral peaks and opens the HOMO-LUMO gap. 



PACS numbers: 72.10.Bg,71.10.-w,73.63.-b,85.30.Mn 

I. INTRODUCTION 

The electron transport through molecular devices has 
gained remarkable interest during last years, primar- 
ily due to experimental adv anc es in creating conduc- 
tive molecule-metal j unctions .1^^ From the experimental 
point of view these systems are very attractive for their 
potential utilization as the next-generation nanometer 
scale building blocks for future integrated circuits ex- 
ceeding up to terahertz operating frequencies. For the- 
orists, the experimental realization of electron transport 
through molecules opens up a new intriguing and chal- 
lenging playground for both theoretical and numerical 
modelling of the underlying physical processes. Under- 
standing these processes at a microscopic level is crucial 
for the future development of molecular electronics. 

Considerable progress has been made to investi- 
gate both steady-stat^2Hi3l ^nd time-dependentPlti^illEll 

transport properties of metal-nanostructure-metal junc- 
tions. As an increasing trend the system is partitioned 
into an explicitely treated interacting region coupled to 
noninteracting electron reservoirs (leads) which act as 
source and sink terminals. However, the partitioning into 
an interacting and a noninteracting part is, in general, 
not well justified due to the long range nature of the 
Coulomb interaction. Recently, there have been some 
advances in calculating transport properties of nanoscale 
junctions while incorporating the electron-elecron inter- 
action in the leads. Perfetto et al.'^^ recently found 
that modelling the electron-electron interaction in low- 



dimensional leads with the Luttinger model the initial 
correlation effects are not washed out in the long-time 
limit and contribute substantially to the steady-state cur- 
rent. Bohr et al¥^ and Borda et al¥^ investigated the 
effects of the lead-molecule interactions in the interact- 
ing resonant level model and showed that it can lead 
to a strong enhancement of the conductance. More re- 
cently these studies hav e been extended to long-range 
lead-molecule interaction^^SEZI, 

Considerable attention has also been devoted to the ef- 
fects of surface polarization (or image charge formation) . 
In Refs. f38ti42l it was shown that polarization effects can 
dramatically change the quasi-particle gap of molecules 
near the metallic surfaces where the dynamical correla- 
tion effects and molecule-lead hopping integrals reduce 
the molecular energy gap across the binding regime from 
gas phase to physisorption. Clearly, this renormalization 
of the molecular levels can have a large impact on the 
transport properties of weakly coupled molecular junc- 
tions. Yet, the question of how the molecule-lead interac- 
tions and, consequently, the formation of an image charge 
affects the ultrafast electron dynamics before a steady- 
state (if any) is reached is still unanswered. The present 
paper wants to address two fundamental issues: what is 
the time-scale to screen molecular charge fluctations in- 
duced by the sudden switch-on of an external bias? And 
what are the scattering processes (or Feynman diagrams) 
relevant for an accurate description of the screening and 
relaxation dynamics? 

To answer these questions we will use the Kadanoff- 
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Baym method which has recently been apphed to both 
finite isolated^^ ''^ and quantum transport systemJ^^^li^ 
and has the merit of preserving all basic conservation 
lawiPSHa. We show that the mean-field Hartree-Fock 
approximation suffers from several limitations in this 
context. Besides being unable to account for dynamical 
polarization effects the Hartree-Fock approximation can 
give rise to "unstable" time-dependent solutions with 
persistent oscillations in density and current. All mean- 
field artifacts disappear when including polarization 
effects in the self-energy, either at the second-Born or 
GW level. These correlated solutions have recently been 
assessed in the Anderson model^*^ and good agreement 
with time-dependent Density Matrix Renormalization 
Group (DMRG) data was found.'^Here we employ them 
for a thorough analysis of the screening versus relaxation 
dynamics as a function of the interaction strength, 
the molecule-lead hopping integrals and the external 
bias. We find that the relaxation time Trei becomes 
shorter when increasing the molecule-lead interactions 
at second-Born and GW level while the screening time 
Tscr is roughly independent on the interaction strength. 
Often, the time-dependent quantum transport simula- 
tions are based on the assumption that Tscr/TVei ^ 1- 
Our results show that the molecule-lead interaction can 
substantially increase this ratio. Another remarkable 
effect of the molecule-lead interaction is that for large 
enough biases the electronic correlations can sharpen 
the spectral peaks and widen the gap between the levels 
of Highest Occupied Molecular Orbital (HOMO) and 
Lowest Unoccupied Molecular Orbital (LUMO). This 
behavior is exactly the opposite of the equilibrium 
behavior and indicates that in the presence of a current 
flow the screening lenghtens the HOMO-LUMO quasi- 
particle life-time and decreases (increases) the ionization 
potential (electron affinity). 

The article is organized as follows. In Section |TI] 
we introduce the model Hamiltonian for quantum 
transport simulations and discuss the exact solution for 
zero molecule-lead hopping integrals. We also give a 
short account of the theoretical background and defer 
the reader to previously published work for details. In 
Section III we analyze the screening versus relaxation 



time and the effect of the formation of an image charge 
in the equilibrium spectral function. Section |IV| deals 
with the short-time dynamics of the lead-molecule-lead 
junction driven out of equilibrium by the sudden switch- 
on of a constant bias while Section |V] deals with the 
long-time dynamics, and in particular with the absence 
of relaxation within HF and the effects of screening in 
the I — V characteristic. The main conclusions are then 
drawn in Section [VII 
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FIG. 1: Image charge model for quantum transport. 

II. IMAGE CHARGE MODEL 

A. Hamiltonian 

To study the image charge effect we c onside r a model 
Hamiltonian that was introduced in Refs.*^^!^. This im- 
age charge model Hamiltonian is displayed schematically 
in Fig. [1] 



H{t)=H,,,,, + H,,,{t) + V ~ ^xN. 



(1) 



The molecular region is modelled by a two-level system 
representing the Highest Occupied Molecular Orbital {H) 
and the Lowest Unoccupied Molecular Orbital (L) with 
energies en and el respectively 



+ UoifiHtnHi + nLtfi-Li) + UHLnnfiL: (2) 



The interaction strengths C/o and Uhl account for the 
intra-level and inter-level electron repulsion. Further- 
more, we used the standard notation hi = X](T=t^ '^i'^ 
for the particle number operator of the molecular level 
i = H,L, where nio- = cl^Ci„ and and Cia are the 
electron creation and annihilation operators. 

The second term in Eq. ([ij describes the left (a — I) 
and right {a = r) leads. 



^ch(i) 



E 

a—l,r — l (T—H 



(3) 



which are modelled as one-dimensional semi-infinite 
tight-binding (TB) chains subject to time-dependent uni- 
form bias voltages W" (t) . The TB parameters hij of the 
chain are chosen so that hij — b for i,j nearest neigh- 
bours and zero otherwise. Finally, c^^^ and Caja are the 
creation and annihilation operators for electrons in lead 
a, site 1 = 1,2,... and spin a. 

The third term in Eq. ([I]) describes the interaction be- 
tween the molecular levels and the TB chains 



^ = E E E A"(cii,c..+cLcai.) 

a=l,r i=H,L <T=t^ 

+ E C^"(«al-l)(iV,„ol-2). 

oc—l,r 



(4) 
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Here A" and [/" are the hopping integrals (proportional 
to the hybridization of the molecular levels) and Coulomb 
interaction strengths between the HOMO/LUMO lev- 
els and the terminal site of lead a. The quantity n^i 
is the particle number operator of site 1 of lead a, 
rial = J2a=^i clio-Caio-j whilc Nraol IS thc total number of 
particle operator of the molecule, iVmoi = fiH + fi^. We 
consider the system initially in equilibrium at zero tem- 
perature, zero bias, = 0, and at half-filling. Then, 
the average density on the lead sites is unity while the 
average density of the HOMO and LUMO levels is 2 and 
respectively. To guarantee the charge neutrality of the 
interacting region we subtracted a positive background 
charge of 1 from hai and of 2 from A^moi- 

This complete the explanation and justification of the 
image charge model (ICM). It can be considered as an ex- 
tension of the interacting resonant level model to study 
molecular excitons and polarization effects. The ICM 
can, of course, be further refined by including interac- 
tions in the leads and a direct lead-lead interaction, and 
can be further generalized to two- or three-dimensional 
leads, more molecular levels, etc. Equation ([!]), however, 
provides the minimal model to study the effects of image 
charges in the non-equilibrium properties of nanoscale 
junctions and in this paper we will not discuss any of the 
aforementioned extensions. 



B. Uncontacted case: Exact solution 

The ICM can be solved exactly for zero hybridization, 
i.e., A*" = a' = 0. In this case the operators fin and 
commute with the Hamiltonian and hence the number of 
electrons on the H and L levels are conserved quantities. 
Let us consider for simplicitv the unperturbed Hamilto- 
nian H obtained from Eq. (hi) by setting the bias W"' to 
zero. All eigenstates of H have the form 



\M,s) = \{e,c]\^s). 



(5) 



Here the -operators create electrons on the molecular 
level j G {H X,H l,L f, ^ 1} and 9j is either equal to 
one or zero depending on what states one likes to occupy. 
The corresponding molecular configuration is specified by 
the collective quantum number M. The state is the 
s-th excited state of the uncontacted leads and has the 
property nj|$s) = 0. For example, a state with two 
electrons in the HOMO-level of the molecule is \H "[,11 1 
,s) = c^.j.c^^|$s). To find the secular equation for the 
|$s) we apply H to \M, s) and find 



H\M, s) 



Hmol + i^ch + ^ f/"(^ial - l)(^mol - 2) 



= [Em + SmA \M, s) 
= EmJM, s). 



\M,s) 
(6) 



where Em is the total energy of the isolated molecule 
with A'liioi electrons satifying the eigenvalue equation 



H^oi\M,s) ^Em\M, s) 



(7) 



while £m,s is the total energy of the uncontacted leads in 
the presence of the potential [/"(A^moi — 2) at the terminal 
sites 



Heh((7)|M,s) 

HchiU ^0)+J2 U"{N^oi - 2)(n„i - 1) 

a 

Sm..s\M, s). 



\M,s) 
(8) 



This potential depends on the strength of the Coulomb 
interaction C/" and on the number N^^i of electrons on 
the molecule. Once we know the electronic configura- 
tion of the molecule, the problem reduces to solving the 
eigenvalue equation Q for a noninteracting TB chain 
with an impurity-like potential at the terminal site. If 
the molecule is charge-neutral, A^moi — 2, this potential 
is zero. However, adding (removing) an electron from 
the charge-neutral molecule gives rise to a potential -I-J7" 
(—[/"). This, in turn, causes a depletion/accumulation 
of charge which is exactly the image charge. 

It is worth stressing that the presence of the lead- 
molecule interaction affects the total energies of the 
charged system, see again Eq. (|8|, and consequently 
changes the addition and removal energies. Consider, 
for instance, the solution for a simple 2-site chain and a 
lead-molecule interaction = U and C/' = (no cou- 
pling to the left lead). It is easy to show that the electron 
afiinity is A = eL+2J7HL+2|&|-2v^(J7/2)2 + while the 
ionization energy is I ^ en + Uq ~ 2\b\ + 2^(C//2)2 6^ 
(see Appendix |A]). The difference A — I reduces with in- 
creasing U and the quasi-particle gap collapses. This 
can also be viewed from another, more general, point of 
view. Consider for simplicity that U" = U for both leads 
and that the intra-molecular interactions Uq and Uhl are 
zero. If the molecule is charge neutral (A^moi = 2) the en- 
ergies of the A'^ and A ± 1 particle ground states (with 
the constraint that the electron is added to or removed 
from the molecule) are given by 



En 
En+1 
En -I 



2eH 



+ fGs(0) 



(9) 
(10) 
(11) 



where we defined £gs{U) to be the ground state energy 
of the Hamiltonian Hch{U) of Eq. Therefore, the 

electron afiinity A and ionization energy / read 



A — En+1 — En 
I = En — En-1 = 



= eL + £gs{U) - £gs{0) 
eH+£Gs{0)-£Gs{-U). 



(12) 
(13) 

Let |$Gs(w)) be the ground state of Hch{U). Then, ac- 
cording to Hellman-Feynman theorenP^I 



d£Gs{u) 
du 



(*Gs(u)|^f^|$Gs(u)), 

du 



(14) 
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FIG. 2: Top left and right panels: The real and imaginary 
part of the dynamical response function. Bottom panel: The 
electron density at the terminal site of the TB chain as a func- 
tion of time when the impurity potential f7 = 0.5 is suddenly 
switched on. The different curves correspond to different val- 
ues of the hopping parameter b = —0.5, —1.0, —1.5, —2.0. 



and therefore 



Sgs{U) - £gs{0) = V / Ki{u) - l]d 

„ "'0 



u. (15) 



From this equation we see clearly how the ground state 
energy depends on the molecular occupation: If we add 
an electron to the molecule we push away charge from the 
first sites of the leads and hence the integral is negative 
and the affinity lowers. On the other hand, if we remove 
an electron from the molecule we attract charge to the 
first sites of the leads and the ionization energy increases. 

The bottom panel of Fig. [2] shows how the image 
charge is built up in the lead. We plot the time-evolution 
of the density at the first site of a semi-infinite chain 
when the impurity-like potential U = 0.5 is suddenly 
switched on at time t — on site 1. The different 
curves correspond to different hopping parameter in the 



lead b = —0.5, —1.0, —1.5, —2.0. By increasing b the 
frequency of the transient oscillations increases and the 
steady state is reached faster. This behavior can be 
easily understood by inspecting the imaginary part of 
the density response function xii(w), top right panel 
of Fig. [2] In Appendix [B] we show that this quantity 
has a maximum at a; ^ 2\b\ which corresponds to 
the oscillation frequency of the density. The width of 
the maximum grows like 2|6| and its inverse gives the 
screening time, i.e., the time-scale for the image charge 
formation. Furthermore, from the top left panel we 
see that Xii(<^ — 0) behaves as 1/5 which is consis- 
tent with the larger induced charge in the long time limit. 



C. Many-body treatment 

The ICM has not exact solution for the contacted case 
and to study it both in and out of equilibrium we use the 
non-equilibrium Green function (NEGF) method based 
on time- propagatio n of the embedded Kadanoff-Baym 
equations.E^E^HlMl xhe basic quantity in the formalism 
is the one-particle Green function 



Gki{z,z') = -i- 



Tr 



{r, 



e */cd2 ^(^)cfc(z)c|(^')]} 



Tr 



(16) 



where we used the notation and c| to denote electron 
annihilation and creation operators either in the molecule 
or in the leads. In the above definition z, z' are the time 
indices on the Keldysh contour c, T is the time-ordering 
operator on the Keldysh contour and Tr{...} signifies the 
trace over the Fock space of all many-body states. The 
Green function G is the solution of the integro-differential 
equation of motion on the Keldysh contour 

[i5^-/i(z)]G(z,z') = '^(^,^')+ y"dz S[G](z,z)G(z,/), 

(17) 

where h{^z") is the Hamiltonian in the one-particle Hilbert 
space, (5(z,z') is the contour delta function and S[G] is 
the self-energy kernel containing all the in formation on 
the many-body and embedding effectJi^E^ Por the pur- 
pose of a practical implementation of the Hamiltonian 
([1]), we divide the system into interacting (C) and non- 
interacting (a) regions and write the single-particle part 
and the interaction part of C as (see Eq. ([l)) 
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FIG. 3: Self-energy diagrams for A' = A'' = 0. A) The only 
nonzero diagrams of the GW self-energy. B) Second order 
self-energy diagrams that are zero. All higher order diagrams 
are also zero. 



Using this notation, the Hamiltonian ([I]) transforms into 



a—l,r ijGot cr—^-i 



(19) 



where C contains the molecular levels and also the 
terminal sites of the leads subjected to the bias voltages 
W^lt). Furthermore, a = L,R are the noninter- 
acting parts of the left and right leads. We choose 
Uo = Uhl = Ulh = 1, M = 0, = -2 and eL = -1 
throughout the rest of this article. 



We will solve Eq.(17l with a Hartree-Fock (HF), 



second-Born (2B) and GW many-body self-energy. The 
quality of the 2B and GW self-energy has recently been 
assessed in the Anderson model^'^ by comparing the time- 
dependent current and density against time-dependent 
DMRG results.!^ Good agreement was found in the pa- 
rameter regime that we discuss below. 

It is instructive and useful for our later analysis to dis- 
cuss the many-body approximations in the uncontacted 
case. Since the number of electrons in the H and L levels 
are conserved quantities in this case, the Green's func- 
tion Gnk = SnkGHH and similarly GLk = SluGll for 
all levels/sites k of the system. The HF approximation 
consists of the first two diagrams in Fig. [s] A). The la- 
bel of the vertices refer to the molecular levels H, L and 
the terminal site of lead a (which is simply denoted by 



1 independently of a). The 2B self-energy is obtained 
by adding to the HF self-energy the first bubble dia- 
gram and the second-order exchange diagram. For zero 
hybridization, however, the second-order exchange dia- 
gram vanishes since it either contains an off-diagonal ele- 
ment of the Green's function (which is zero) or a product 
of lesser and greater diagonal element of the molecular 
Green's functions, e.g., G^^^G^^. Having the equilib- 
rium system two electrons in H and zero in L it must be 
Gf^ = ^Hff = 0. Consequently, only the first bubble 
diagram survives and the 2B self-energy diagrams are all 
displayed in Fig. [sjA). Note that if the external vertices 
of the bubble self-energy diagram lie on the terminal site 
of the leads [second diagram of Fig. [3] B)] the diagram 
vanishes. This is a direct consequence of the fact that 
the polarization diagram with or L as vertices is pro- 
portional to the product G^^G^^ or G^^G^^ which is 
zero. For the same reason the first diagram of Fig. |3] B) 
is also zero. The physical origin of this result is that one 
cannot create particle-hole excitations on the molecules 
without changing the number of electrons in the H/L 
levels. The many-body self-energy in the GW approxi- 
mation is S = iGW where the screened interaction W is 
approximated as a geometric series of bare polarization 
diagrams connected by interaction lines. Since the only 
bare polarization diagram is the particle-hole propagator 
going from 1 to 1, the GW approximation coincides with 
the 2B approximation. In our ICM there is no direct in- 
teraction between two electrons on the terminal site of 
the leads. We therefore expect the 2B and GW approxi- 
mations to perform similarly for small hybridizations. 



III. COMPETING TIME-SCALES AND 
SPECTRAL PROPERTIES 

In the previous Section we have seen that there is a 
characteristic screening time to build up charge after 
the addition or removal of an electron to or from the 
molecule. In the case that the molecule is contacted to 
the leads there is another time-scale that plays a role. 
This is the relaxation time to disperse the excess charge 
on the molecule into the leads. It is the ratio between 
the screening time and the relaxation time that tells us 
how the system behaves under non-equilibrium condi- 
tions. The aim of this Section is to extract these time- 
scales from the equilibrium properties of the contacted 
system and to analyze the effects of the screening on the 
equilibrium spectral function. This will help us to gain 
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FIG. 4: 



Green function G^j^ {t, 0) for a) 



HF and 2B with 



b = -0.6, [/ = 1 and A = 0, b) HF and 2B with b = -O.i 
U = 1.0 and A = -0.2 
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FIG. 5: Green function Gl^{t,0) for HF and 2B with b = 
— 1.0 (panels a) and b) ) and b — —0.6 (panels c) and d)) and 
A — —0.2 and different values of U. 



insight in the more complicated case of quantum trans- 
port discussed in the next Section. The analysis will be 
carried on using many-body Green's function methods 
since the contacted case is no longer analytically solv- 
able. 



A. Screening and relaxation times 

In this Section we study the response of the system 
to the sudden addition or removal of an electron on the 
molecule within the HF, 2B and GW approximations. 
The response of an added electron is encoded into the 
and Green functions which we can calculate within 
these many-body approximations both in real time and 
in frequency space. For instance, the LUMO Green's 
function Gl^{t,0) = -i{cLit)c{{0)) gives the probabil- 
ity amplitude of finding a particle on the LUMO level 
at time t after being created at time 0. In Fig. |4] we 
plot the real part of this quantity. This Green function 
oscillates with a characteristic frequency equal to the ad- 
dition energy of an electron to the LUMO level. In the 
panel a) of Fig. |4] we compare the HF and 2B results 
for A' = A'' = 0, f/' = 0, f/'' = C/ = 1 and b = -0.6. 
The correlated 2B curve exhibits a short transient with 
a characteristic time-scale ^ 1/5. This transient has to 
be attributed to the build up of the image charge and its 
duration is consistent with the previous analysis of Fig. 
[2j Note that no transient is visible in the HF approxima- 
tion which therefore fails to describe the formation of the 
image charge. The interplay between the screening time 
and the relaxation time can be investigated by contacting 
the molecule to the leads. In the panel b) we consider the 
same parameters as in panel a) except for A'' = A = —0.2, 
i.e. the molecule is contacted to the right lead. The main 
difference to the previous case is that both the HF and 
2B curves are damped (relaxation). Similarly to the un- 



contacted case, there is no evidence of screening in the 
HF approximation. 

Let us now address in more detail the dependence of 
the relaxation time on the molecule-lead interaction U. 
In most physical situations the band-width of the metal- 
lic leads is much larger than the molecule-lead coupling, 
b ^ X. Then, for small values of U the relaxation time 
is proportional to t^ci ^ T^^ ~ l^l/''^^ • This time- 
scale depends both on the molecule-lead coupling and 
the lead hopping. On the other hand, the screening time 
Tscr ^ l/l^l is a property of the lead only and the ratio 
Tj-ci/Tgcr ^ {b/X)^ is always larger than unity. If the in- 
teraction U becomes comparable to or larger than b then 
this analysis is not valid anymore. In the Appendix [C] 
we show that already at the HF level U renormalizes the 
relaxation time according to Tj-ei ~ r~^(l — CUY where 
C is a positive real constant weakly dependent on U for 
small values of U . This is illustrated in the top panels 
of Fig.([5| where we display the real part of G^j^{t, 0) for 
b = -LO and A'' = A = -0.2 at HF and 2B level. We 
see in both cases that by increasing the molecule-lead 
interaction U we lower the relaxation time. The renor- 
malization of the lead coupling (or the embedding SE) 
also affects the positions of the molecular quasiparticle 
levels. The renormalization leads to a small upward shift 
of the LUMO level and a downward shift of the HOMO 
level, i.e., a slight opening of the HOMO-LUMO gap. 
This is clearly visible in Fig. [5^, where we see a slight 
increase in the frequency of the LUMO oscillation when 
we increase U . In panel b) for 2B, on the other hand, 
we see a much more drastic decrease of the oscillation 
frequency due to the image charge effect which HF fails 
to describe properly. 

We note that the small upward shift within the HF 
approximation of the LUMO level with increasing U can 
lead to an increase of the relaxation time when the level is 
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close to the band edge. This is because the upward shift 
pushes the LUMO level close to the band edge where the 
imaginary part of the embedding self-energy decreases 
rapidly and compensates the renormalization introduced 
by the interaction U . In this case, the spectral peak de- 
scribing the position and lifetime of the molecular quasi- 
particle level becomes also highly asymmetric and non- 
Lorenzian which leads to a non-exponential decay of the 
Green function in real time. These features are illus- 
trated in Fig. [5]; where we consider the case of lead 
hopping b = -0.6. The LUMO level for t/ = is located 
at 1 which is close to the band edge of 2\h\ = 1.2. An 
increase of U to 1.3 pushes the level very close to the 
band edge and we then see a corresponding increase in 
relaxation time with a non-exponential decay. In the case 
of 2B (panel d)) the image charge effect pushes the level 
inwards, away from the band edge, and we see that the 
relaxation time again decreases with increasing U . Com- 
paring panel d) to panel b) we see further that increasing 
the lead hopping h leads to a slight decrease of the image- 
charge effect (frequency change for U ^ Q) and increase 
of the r elaxa tion time, in agreement with the analysis of 
Section II B and the relation t^c\ ^ ^ l^l/'^^- 



The difference between the HF and the correlated re- 
sults in real time translates into a different spectral struc- 
ture in frequency space. In Fig. [6] we display the quasipar- 
ticle spectral functions (see Eq. 20 ) for the LUMO level. 



All{lo), corresponding to the Green functions G^^(t,0). 
This is done for the 2B approximation using h — —0.6 and 
different values of U . For the uncontacted case the 2B re- 
sult coincides with the GW result, as discussed in Section 
II C[ The corresponding spectral function for U = 1.0 is 
displayed in the left panel while in the right panel we 
have A = —0.2 and we plot the spectral function for dif- 
ferent values of U . The very fast oscillations in the left 
panel are due to the finite time interval in the Fourier 
transform. They are not present in the right panel due 
to damping of the Green function in the contacted case. 
Besides the main peak located at the electron affinity 
we observe a shoulder of width A « 4|6| at higher en- 
ergies. At finite hybridization A = —0.2 this shoulder is 
smoother and partially merges with the main peak. The 
shoulder originates from the particle-hole continuum of 
excitations induced by the sudden addition of an electron 
to the LUMO state. They are these excitations which al- 
low for the dynamical screening of the extra charge on the 
molecule. In mathematical terms the shoulder arises by 
Fourier transforming the initial transient of the 2B curve 
in Fig. [4]and[5j Since no transient was observed in HF, 
the HF spectral function will consist only of a main peak 
at the electron affinity. Both 2B and GW incorporate 
the correct physics through the polarization diagram of 
Fig. [Sj'V, which nicely illustrate how an extra electron on 
the LUMO can excite a particle-hole on the terminal site 
of the leads. For small hybridizations the polarization is 
approximatively equal to the response function of Fig. [2] 
which explains the width 4|6| of the shoulder. We further 
see in the right panel of Fig. |6]that while the peak moves 
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FIG. 6: LUMO spectral function in the 2B approximation for 
h = -0.6, U = 1.0 and A = (left panel), h = -0.6, A = -0.2 
and (7 = 0, 1.0, 1.3 (right panel). 



leftward with increasing U the width of the plateau re- 
mains roughly constant at 4|6|. The screening time is 
therefore independent of the molecule-lead interaction. 

The message to take home is that the molecule-lead 
interaction have a large impact on the ratio Tici/tsci- and, 
in principle, can turn it to be smaller than one. This 
kind of exotic situations would occur in leads with flat 
bands as, e.g., those modeled by Tasaki.^'* In most metal- 
lic systems this is not the case and in the remainder of 
this paper we will study the regime Troi/Tgcr > 1. 



B. Equilibrium spectral function 

In this Section we investigate the effects of screening 
on the spectral features of the molecule in equilibrium. 
We calculate the molecular spectral function A^o\{ijj) as 
a sum of the projected spectral components Aii{Lo) as 



(20) 



i=H,L 



for J7' = a' = and for zero and finite hybridization 
A*" = A with the right lead. The results are shown in Fig. 
[7] for A = and in Fig. |8]for A = -0.2. 

Let us start by analyzing the performance of the HF 
approximation. The first observation is that for A = 
the HOMO-LUMO gap and the intensities of the peaks 
remain unchanged as the interaction strength — U 
increases. This can easily be understood from the explicit 
form of the HF HOMO and LUMO energies 



HF 



{en 



U) 
U) 



uhUo 
ulUo - 



- 2nLUHL 



2niU, (21) 
2niU. (22) 



At half-filling the average density ni — nir of the right 
terminal site is 1/2 and hence the dependence on U can- 
cels off. In the case of finite hybridization A = —0.2 (Fig. 
|8]), the HF peaks shift slightly outwards and broaden 
due to the renormalization of the embedding self-energy 
(or equivalently, the renormalization of the hybridization 
A — > A + C/G^]^, see Appendix[c|. It is then clear that for 
A = the intensities do not change since Ghi = 0. Simi- 
lar renormalization effects has been observed in Ref. ESI 
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FIG. 7: HF and GW equilibrium spectral functions Amoil'^) 
for A = 0, 6 = -1.0 and U = 0, 0.5, 1.0, 1.5. The vertical lines 
indicate the exact and 2B(HF) peak positions. 
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FIG. 8: HF, 2B and GW equilibrium spectral functions 
Amoi{(^) for A = -0.2, b = -1.0 and [7 = 0, 0.5, 1.0, 1.5. 
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FIG. 9: The real and imaginary components of the 2B(HF) 
self-energy for the molecular HOMO level with U — 0.5 (left 
panel) and ?7 = 1.5 (right panel). The rest of the parameters 
are Uhl — Uo = — b ~ 1 and A = 0. 



the charged system can lower its energy by exciting par- 
ticles from occupied to unoccupied levels of the charge- 
neutral system. The resulting effect is to accumulate or 
deplete charge in the neighborhood of the terminal site, 
i.e., to screen the excess charge of the molecule. Note also 
that in the case A 7^ 0, the coupling to the particle-hole 
continuum provides an extra channel for quasi-particle 
scattering and induces quasiparticle broadening to the 
spectral peaks. The differences between the uncontacted 
and contacted spectral functions must be attributed to 
charge transfer processes and the consequent formation 
of image charges in the molecule. This molecular polar- 
ization effect was recently found to reduce the HOMO- 
LUMO gap even further.ESi 

To assess the quality of the correlated approxima- 
tions and the importance of self-consistency we display 
in Fig. [7] the position of the exact H/L peak (calculated 
from the Hellman-Feynman theorem) as well as the po- 
sition of the peaks as obtained from a one-shot 2B cal- 
culation with HF Green function ( denoted with 2B(HF) 
). As can be seen from Fig. [7] the GW results arc m very 
good agreement with the exact ones. The position of 
the spectral peaks in the correlated approximations are 
obtained from the quasiparticle equation 



,HF 



Re{l]«(a;)} =0, 



(23) 



In the HF approximation the self-energy of the H/L lev- 
els couples only to the density at the terminal site of the 
lead and thus misses entirely the particle-hole coupling 
responsible for the screening. 

The situation is radically different in the correlated 
2B and GW approximations. In both cases the HOMO- 
LUMO gap, corresponding to the difference A — I be- 
tween the electron affinity and the ionization poten- 
tial, narrows in agreement with the discussion of Section 
II B The added/removed electron and its image charge 



bind together, thereby decreasing/increasing the addi- 
tion/removal energy. The stronger is the interaction U 
and the larger is the gap reduction. In the 2B and GW 
approximations the added/removed electron couples not 
only to the density but also to the particle-hole contin- 
uum of the lead. It is through this latter coupling that 



where S^(cli) is the retarded many-body self-energy pro- 
jected onto the i = H,L molecular level. In Fig. |9] we 
display the graphical solution of Eq. (23) with 2B(HF) 
self-energy and i = H. For this plot we have chosen 
Uhl = Uo — ~ 1 and U — 0.5 (left panel) and 
U = 1.5 (right panel). Already one iteration of the self- 
consistency cycle captures the correct trend. The zero of 
Eq. ( 23 1 moves toward higher energies with increasing U. 



An analogous calculation for the LUMO level shows that 
the zero moves toward lower energy. In conclusion, the 
inclusion of polarization effects into the self-energy has 
two main effects in equilibrium: the redistribution of the 
spectral weight due to particle-hole excitations (satellite 
spectrum) and the collapse of the HOMO-LUMO gap. 
As we shall see, the situation is radically different out of 
equilibrium. 
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FIG. 10: a) Time-dependent right current for U = 
0.0,0.5, 1.0. b) Ground state (GS) and steady state spectral 
function of Eq. (20 1, c) Time-dependent densities nir(t) and 

d) Time-dependent HOMO and 
In all the plots the sim- 



nii{t) at the terminal sites. 
LUMO densities nnit) and nL(t). 
ulations have been performed within the HF approximation 
with bias = -W = 0.8. 



FIG. 11: a) Time-dependent right current for U = 
0.0,0.5,1.0. b) Ground state (GS) and steady state spectral 
function of Eq. (201. c) Time-dependent densities n-irit) and 

d) Time-dependent HOMO and 
In all the plots the sim- 



nii{t) at the terminal sites. 
LUMO densities nnit) and nL{t). 
ulations have been performed within the HF approximation 
with bias = -W = 1.2. 



IV. QUANTUM TRANSPORT: SHORT-TIME 
DYNAMICS 

In order to investigate the short-time transport prop- 
erties of the system of Fig. [l] we consider A' = A'' = A 
and [/' = = U . We will analyze the transient dynam- 
ics after the sudden switch-on of a bias — —W^ = W 
in the leads. Note from Eqs. (18 1 and (19) that the 



bias is applied also to the terminal (interacting) sites 
of the leads. We will refer to the left /right current as 
the current flowing through the left/right interacting- 
noninteracting interfaces correspondingly. In all simu- 
lations we set A = —0.2 and hence work in the weak 
tunneling regime to highlight correlation effects. 



A. HF approximation 

In Figs. [To] and [TT] wc show the time-dependent cur- 
rents (panel a), ground state and nonequilibrium steady 
state spectral fimctions (panel b), terminal site densi- 
ties (panel c) and HOMO/LUMO densities (panel d) 
for the HF approximation with molecule-lead interaction 
U = 0.0,0.5, 1.0. In Fig. [10] we consider the "small" bias 
case W^' = —W^ = 0.8 for which the equilibrium H/L 
levels e^^L ~ "pl remain outside the bias window while 
the bias is set to W'' 



in Fig. 



11 



-W = 1.2 so that the 
equilibrium H/L levels lie inside the bias window. 

For zero molecule-lead interaction, U — 0, and small 
bias the current flowing through the system is almost 
zero, see Fig. [TO^ , in agreement with the fact that the 
H/L levels are outside the bias window. A finite current 
instead sets in for large bias, see Fig. [TT^ . The physics 



is here very similar to that of the non-interacting reso- 
nant transport regime. On the other hand, the current 
increases substantially at finite U for small bias. Fur- 
thermore, increasing U the frequency and the amplitude 
of the oscillations in the current and density becomes 
larger. We recall that in the HF approximation the equi- 
librium quantities are fairly independent of U. These re- 
sults show that at finite bias the situation is completely 
different. 

To understand the differences between the equilibrium 
and the non-equilibrium case we observe that the gap 
in the non-equilibrium spectral function reduces consid- 
erably at finite U. For instance Fig. [T0| 3 shows that for 
U = 0.5 and U — 1.0 the H/L levels have already entered 
the bias window [—0.8,0.8]. To trace back the physical 
origin of this effect we write the HF energies of the H/ L 
levels 



= en -2U + UanH + 2UHLnL + '2U[nir 



nil] 



= eL-2U + UonL + 2UHLnH + 2U[nu.+nu], 

where we took into account that the molecule is now con- 
nected to both leads. The terms containing an explicit 
dependence on U cancel off since the sum of the termi- 
nal site densities, nir{t) -(- ni;(t), remains roughly at its 
ground state value during the entire time evolution, see 
panels c). Thus, it is not the lead polarization which af- 
fects the level positions but rather the polarization of the 
molecular region, i.e., the difference rt// — n^. The panels 
d) indicate that the molecular polarization increases as 
U becomes large. This analysis shows that in the HF ap- 
proximation the reduction of the gap induced by U has 
the same nature observed earlier^ and has nothing to 
do with the image charge effect. However, as we will see 
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FIG. 12: a) Fourier transform of the transient current of Fig. 
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with U = 1.0. b) Ground state and steady-state spectral 
function j4moi(i^) for U — 1.0. c) Time-dependent density 
matrix components G'L,i;(t, t^) and Gu,L{t,t^) for U = 0.0 
and U = 1.0 d) HF time-dependent density matrix compo- 
nents GH,ii{t,t+) and Gii,Hit,t+) for U = 0.0 and U = 1.0. 
In all the plots the simulations have been performed within 
the HF approximation with bias = —W^ = 1.2. 



later, this effect already has a big impact on the resulting 
current-voltage characteristics. 

The main frequency of the oscillations in the transient 
density and current originate from the electronic transi- 
tions from the left electrochemical potential /x' = ep + W'' 
to the LUMO level and also from the HOMO level to the 
right electrochemical potential /i'' = ep + (for the 
symmetric bias considered here these transitions have 
the same energy). This can easily be verified by cal- 
culating the discrete Fourier transform of the transient 
current, /(w). In Fig. 12 we show /(cj) for U ~ 1.0 and 
the large bias case W^~= — 1.2 (panel a) along 

with the ground state and nonequilibrium steady-state 
spectral function (panel b). The Fourier transform I{uj) 
exhibits a sharp peak at w ~ 1.0 with a smearing to- 
wards lower frequencies down to uj w 0.2. The smearing 
is a direct consequence of including the transient part 



of I{t) in the Fourier transform. The value of e|^^^ is 
=F1.0 in equilibrium while it is about ^0.2 at the steady 



state, see Fig. [12)3. As the HOMO-LUMO gap collapses, 
the transition energy between the left /right electrochem- 
ical potential and the LUMO/HOMO level changes from 
0.2 to 1.0. The aforementioned smearing towards low 
frequency is the fingerprint of the dynamical renormal- 
ization of the transition frequency. Another consequence 
of the collapse of the steady-state gap with increasing 
U is that e^^^ moves further away from /i'/'' where the 
density of states has a square-root divergence (resonance 
condition). This is clearly illustrated in Fig. [TI| 3. The 
further away the levels are from resonance the harder it 
is for electrons to tunnel, which in turn implies a larger 
oscillation amplitude and a smaller average current. 
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FIG. 13: Time-dependent density in the non-interacting part 
of the right lead for U = 0.0 (top panel) and U = 1.0 (bottom 
panel). The simulations have been performed within the HF 
approximation with bias — —W^ = 1.2. Site number 2 
corresponds to the first noninteracting site in the right lead. 



The transient oscillations are also visible in the off- 
diagonal components of the time-dependent density ma- 
trix, Gij{t) = Gij{t,t^), which is displayed in Figs. 12; 
and [1^ for [/ = and f7 = 1.0. The component 
GL,ii{t) and Gii^L{t) oscillate with the same main fre- 
quency as the current and densities. The same holds 
true for GH,ir{t) and Gir^nit) (not shown). On the con- 



trary GH,ii{t) and Gu^nit) have a very weak high fre- 
quency component superimposed to the main frequency. 
Initially the HOMO level is fully occupied and electronic 
transitions from /i' to are blocked. Similarly, the 
LUMO level is initially empty, so there are no electronic 
transitions from fi^ to e|[^. As the time passes, however, 
the HOMO occupation decreases while the LUMO occu- 
pation increases and these transitions become possible. 
They are located around 1.4 and 2.0 and can be seen in 
the Fourier transform of the current (the current is indeed 
given in terms of off-diagonal elements of the density ma- 
trix). Even though present, the transitions between the 
HOMO level and the LUMO level are extremely small 
since there is no direct hopping between the two levels. 

The sudden switch-on of the bias gives rise to den- 
sity shock waves in the leads with features similar to the 
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FIG. 14: a) Time-dependent right current for U = 
0.0,0.5, 1.0. b) Ground state (GS) and steady state spectral 
function of Eq. (20 1, c) Time-dependent densities nir(t) and 

d) Time-dependent HOMO and 
In all the plots the simula- 



nt (t) at the terminal sites. 
LUMO densities nnit) and nL{t). 
tions have been performed within the 2B approximation with 
bias W' = -W = 1.2. 



FIG. 15: a) Time-dependent right current for U = 
0.0,0.5,1.0. b) Ground state (GS) and steady state spectral 
function of Eq. (201. c) Time-dependent densities nir{t) and 

d) Time-dependent HOMO and 
In all the plots the sim- 



nii{t) at the terminal sites. 
LUMO densities nnit) and nL{t). 
ulations have been performed within the GW approximation 
with bias = -Vy = 1.2. 



density at the terminal sites. In Fig. [13] we show the 
transient dynamics of the HF density in the noninter- 
acting part of the right lead for U = 0.0 (top panel) 
and for U — 1.0 (bottom panel) when the bias voltage is 
— —W^ — 1.2. The shock wave reaches site j after a 
time j /vp where in our case the Fermi velocity vp — 2b. 
No matter how far site j is the density at this site ex- 
hibits damped oscillations whose initial amplitude and 
relaxation time is independent of j and increases with U. 



imation attenuates the gap opening compared to the 2B 
approximation, but the sharpening of the peaks is well 
visible also in this case. The gap opening in the out-of- 
equilibrium system has never been reported before and, 
as we shall see below, has profound consequences on the 
I ~ V curve. 



QUANTUM TRANSPORT: LONG TIME 
DYNAMICS 



B. Correlated approximations 

The inclusion of correlations changes considerably the 
physical picture. Let us focus on the large bias case 
= —W^ = 1.2 and calculate the same quantities as 
in Fig. [TT] but within the 2B and GW approximation. 
The results are displayed in Fig. 14 and 15 respectively. 
The first important feature is that the relaxation time is 
much shorter than in the HF case due to the many-body 
broadening of the HOMO and LUMO levels, see panels 
b). In the same panels we also show the ground state 
(GS) spectral function for the same values of U. As ex- 
pected the GS gap between the HOMO and LUMO peaks 
reduces with increasing U due to the image charge effect. 
In the bias ed system for ?7 = the bias dependent gap 
closinj^l^ brings the levels so close to each other that 
we can observe only one very broad peak. Interestingly 
and surprisingly, the effect of increasing U in the biased 
system is to open the gap and to sharpen the spectral 
peaks. In the 2B approximation with molecule-lead in- 
teraction U = 1.0 the nonequilibrium steady-state gap is 
even larger than the ground-state gap. The GW approx- 



In this Section we investigate the effects of the image 
charge on the long-time dynamics of the lead-molecule- 
lead system within the HF, 2B and GW approximation. 
As we shall see a non-trivial post-transient dynamics de- 
velops at the HF level. The inclusion of correlations does 
always bring the system in a steady-state regime. We 
will show how this regime is attained and calculate cur- 
rent and densities in the steady state for different bias 
voltages and molecule-lead interaction. 



A. HF approximation and post-transient dynamics 

We focus on the large bias regime and strong molecule- 
lead interaction U — 1.0. In the previous Section we 
showed that current and densities seem to relax after the 
transient behavior induced by the sudden switch-on of 
a bias voltage. However, extending further the propaga- 
tion time-window something unexpected occurs. We find 
that the steady state is metastable and oscillations with 
increasing amplitude develop to then stabilize in a peri- 



odic state. In Fig. 16 we display long-time simulations 
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FIG. 16: Time dependent right current (top panel) and total 
number of particles in the molecule (bottom panel) for bias 
voltages in the range [0.7, 1.2] and molecule-lead interaction 
U = 1.0. 



of the right current I{t) (top panel) and the total num- 
ber of particles in the molecule N^oi{t) (bottom panel) 
for bias voltages in the range [0.7, 1.2]. In this range the 
equilibrium HOMO and LUMO levels lie in the bias win- 
dow. The frequency of the oscillations increases as the 
bias voltage is increased, which is a clear indication that 
the dominant transitions are those between the leads and 
the molecular levels. 

In Fig. [it] we display the time-dependent left and 
right currents (panel a) as well as the terminal-site den- 
sities (panels c) and molecular densities (panel d) for 
= —W^ = 0.8. According to these results the 
post-transient periodic state corresponds to a sequence of 
charge blockades with opposite sign of the electron-liquid 
acceleration (time-derivative of the current) between two 
consecutive blockades. The oscillations are therefore due 
to a charge sloshing between the molecular levels and 
the terminal sites. The metastability of a steady-state 
solution in which current and densities are given by the 
average value of the time-dependent results is due to the 



FIG. 17: a) Time-dependent left, right and total current, b) 
Fourier transform of the total current, c) Time-dependent ter- 
minal site densities, d) Time-dependent HOMO and LUMO 
densities and the total number of particles in the molecule. In 
all panels U = 1.0 and the bias voltage is W'' = —W = 0.8. 

combination of the constant flow of electrons from left to 
right and the self-consistent nature of the Hartree-Fock 
potential. Finally we emphasized that the amplitude of 
the ac current superimposed to the dc current depends 
on where the current is measured. 

In Fig. [TTj j we report the Fourier transfrom of the total 
current, /tot(w). The main peak at a; « 0.4 is smeared 
out toward higher frequencies up to w « 0.6, indicat- 
ing the occurrence of electronic transitions between lev- 
els whose position changes dynamically in time. We also 
observe higher frequency satellites arond uj « 1.2, 1.8. 
These satellites occur exactly at the positions of odd har- 
monics of the main frequency. The absence of even har- 
monics is due to the fact that the external driving field 
is an odd function in space. 

The persistent oscillatory behaviour reported in this 
Section is most likely an artifact of the HF approxima- 
tion and, as we shall see in the next Section, disappears 
in the 2B and GW approximations. Within HF the sys- 
tem knows only the instantaneous density and there is no 
damping mechanism to wash out the oscillations. These 
oscillations are sustained by the finite bias voltage and 
originate from the instantaneous Coulombic feedback. 

B. Steady-state properties: HF, 2B and GW 
approximation 

In Fig. [18] we show the HF time-dependent currents 
and the resulting I — V characteristic (bottom right 
panel) for different interaction strengths U — 0.0, 0.5, 1.0. 
Since the HF currents for f7 = 1.0 do not attain a steady 
state for large enough bias, the I — V characteristic is 
in this case calculated with the dc part of the current 
(average value). The inclusion of the molecule-lead in- 
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FIG. 18: HF time-dependent right current for different in- 
teraction strengths U = 0, 0.5 (top left and right panels) and 
U = 1.0 (bottom left panel). The I — V curves extracted from 
the long-time limit are displayed in the bottom right panel. 



FIG. 20: Ground state (GS) and steady-state spectral func- 
tions in the HF approximation for U = 1.0 and bias W'' = 
— 0.55. a) Full spectral function of the interacting re- 
gion, b) Spectral function on the terminal site of the left 
lead, c) HOMO and LUMO spectral functions, d) Spectral 
function on the terminal site of the right lead. 




0.2 _°NJ% 



FIG. 19: Time-dependent number of electrons on the 
molecule, N^oi{t), versus the applied bias voltage. HF ap- 
proximation with U — 1.0. 



teraction deforms the I — V characteristics dramatically. 
Firstly, increasing the interaction strength, the threshold 
is shifted towards smaller bias values. Secondly, increas- 
ing the interaction strength up to [/ = 1.0 gives rise to 
an extra step in the I—V curve. The shift of the I — V 
step towards smaller biases is related to the gap closing 
mechanism which in the HF approximation is entirely 
due to the intramolecular interactions Uq and Uhl-: see 
Section |IV3 

The extra step in the HF I — V curve (bottom panel of 



Fig. 18 1 corresponds to a charged state of the molecule. 



In Fig. 19 we plot the number of particles (per spin) 
in the molecule, N^^i, for interaction U = 1.0. There 
exists a narrow window of applied biases — —W^ £ 
[0.55,0.6] for which Nj^oi ~ 1-35. We have also checked 
(not shown here) that this window can be extended by in- 
creasing the molecule-lead coupling A. The excess molec- 
ular charge produces a Hartree barrier on the terminal 



sites which prevents the current to increase, see plateau 
in the I—V curve for = 1. As the bias becomes larger 
electrons gain enough energy to overcome the barrier and 
the current increases again. 

The excess charge on the molecule changes also the 
spectral function. In Fig. 20 we plot the full spectral 



function of the interacting region as well as the local 
spectral functions of the HOMO, LUMO and the termi- 
nal sites in the ground and steady state for U = 1.0 and 
for bias VF' = —W^ = 0.55 within the HF approxima- 
tion. The HF spectral function of the charged molecule 
exhibits two sharp structures close to the left and right 
band edges (they are separated by — — 1.1). The 
induced Hartree barrier pushes electrons away from the 
terminal sites and gives rise to well localized hole states. 
The structure of the peaks is indeed similar to that of 
a split-off state (anti-bound state) which forms in the 
presence of an external positive potential at the end- 
site of a semi- infinite chain, see Appendix [Xj In our case 
this potential is v" — VF" -I- wh with Hartree potential 
vu = 2i7(iV„,oi - 1) « 0.7. 

The formation of the additional step in the I—V curve 
is probably another artifact of the HF approximation. In 
Figs. [21] and [22] we show the transient and steady-state 
currents for U = 0,0.5 and U = 1.0 and bias voltage in 
the range [0, 1.2] within the 2B and GW approximations. 
Like for the HF approximation the onset of the current 
is shifted towards smaller bias values when U increases. 
However, this effect is more pronounced in the 2B and 
GW approximations which properly incorporate dynam- 
ical polarization effects to account for the formation of 
the image charge. Another effect of correlations is to 
smoothen the onset, in agreement with the appearence 
of a particle-hole shoulder in the spectral function, see 
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FIG. 21: 2B time-dependent right current for different in- 
teraction strengths U = 0, 0.5 (top left and right panels) and 
U = 1.0 (bottom left panel). The I — V curves extracted from 
the long-time limit are displayed in the bottom right panel. 
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FIG. 22: GW time-dependent right current for different in- 
teraction strengths U = 0, 0.5 (top left and right panels) and 
i7 = 1.0 (bottom left panel). The I — V curves extracted from 
the long-time limit are displayed in the bottom right panel. 



Section HnXl 

It is important to disentangle the scattering-induced 
brodeaning due to the intramolecular interactions f/g and 
Uhl from that due to the molecule-lead interaction U. 



In Fig. 23 we plot the full spectral functions of the inter- 
acting region within the HF, 2B and GW approximations 
for two different values of J7 = and U = 1.0. For U = 
and within 2B and GW there is a consistent broaden- 
ing of the HOMO/LUMO peaks when these levels enter 
the bias window. This effect was reported previously in 
Refs. manddll However, for f7 = 1 the 2B and GW spec- 
tral functions do not get broader as they enter the bias 
window. The peaks preserve their shape and the HOMO- 
LUMO gap starts to open up. The molecule-lead interac- 
tion has an effect opposite to that of the intramolecular 
interaction on the broadening and the many-body shift 




FIG. 23: HF, 2B and GW spectral functions for different bias 
voltage and lead-molecule interaction U = 0.0 and U = 1.0. 



of the spectral peaks. This is a very important result 
according to which image charge effects in the biased 
system contribute to lenghtcn the HOMO/LUMO quasi- 
particle lifetimes and decrease (increase) the ionization 
potential (electron affinity). 



VI. CONCLUSIONS 

In conclusion, we provided a thorough analysis of 
the effects of the dynamical formation of image charges 
at the interfaces between a molecule and the metallic 
leads under non-equilibrium conditions. The analysis has 
been carried out within the embedded Kadanoff-Baym 
method using fully self-consistent many-body approxi- 
mations at the HF, 2B and GW level. The mean field 
HF approximation fails to capture the polarization effects 
both in and out of equilibrium. As a consequence, the 
equilibrium molecular levels are not renormalized while 
out of equilibrium the renormalization is solely due to 
the intramolecular interactions. We pointed out that the 
shortcomings of the HF approximation are also at the 
origin of other unphysical effects. There exists a finite 
range of applied biases for which the molecule is artifi- 
cially charged. This causes a depletion of the electron 
density at the interfaces and prevent the current to in- 
crease as the bias becomes larger (plateau in the I—V 
characteristic). Furthermore, for large enough bias and 
molecule-lead interaction the molecular system does not 
relax in the long time limit. We reported the occurrence 
of the undamped oscillations in current and densities. 
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These oscillations correspond to a charge sloshing be- 
tween the molecular levels and the terminal sites. 

To cure the problems of the mean-field theory we re- 
sorted to the 2B and GW approximations. In both ap- 
proximations the self-energy contains polarization dia- 
grams which correctly account for the screening of the 
charged molecule and hence are suited to describe the 
formation of image charges. In all situations considered 
we did not observe a plateau in the I — V characteris- 
tic nor the absence of relaxation. An important finding 
of our analysis is that by increasing the molecule-lead 
interaction the ratio between screening time and the re- 
laxation time changes and the screening time is primarily 
determined by the properties of the lead. As expected, 
the 2B and GW equilibrium HOMO-LUMO gap closes 
when increasing the molecule-lead interaction. Thus, 
the onset of the current in the I ~ V characteristic is 
shifted to lower biases as compared to a non-interacting 
or mean-field calculation. Another remarkable effect per- 
tains the molecule spectral properties as a function of the 
applied bias. In equilibrium the molecule-lead interaction 
is responsible for the reduction of the HOMO-LUMO 
gap and for a substantial redistribution of the spectral 
weight to the satellites induced by the electron correla- 
tions. Increasing the bias the situation changes. For zero 
molecule-lead interaction the HOMO and LUMO peaks 
near each other and considerably broaden when they en- 
ter the bias window. The effect of the molecule-lead in- 
teraction is to keep the spectral peaks sharp and to open 
the HOMO-LUMO gap. This effect is therefore exactly 
the opposite of that generated by the intramolecular in- 
teractions. All this phenomenology clearly shows the im- 
portance of a proper description of electron correlations 
in time-dependent and steady-state quantum transport. 
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FIG. 24: Electronic configuration for a two-site lead. For 
an extra electron on the molecule there are three energy 
eigenspaces for the lead. Note that the parallel-spin electron 
states do not contribute in the response properties since the 
bias preserves the spin orientation. 



M ~ GS^ the molecular configuration with the extra 
electron. At half-filling the right lead has two electrons 
and the eigenstates |GS^s) of Eq. (jsj) are displayed in 
Fig. [24] Their energy is Eq = 2si, £i = ei + £2 and 
£2 = 2e2. In a similar manner we can calculate the 
lead eigenenergies for the molecule with an electron less. 
The resulting ionization potential and electron affinity 
are I ^ - E3 = en + Uq - 2\b\ + 2^{U/2y+b^ and 
A E5_- E4 = tL^+ 2Uhl + 2\b\ - 2VW2F+62, see 
Eqs 



2Uhl + 2\b\ 

(12) and (13 1. These energies correspond to the 



renormalized energies of the HOMO and LUMO level. 
We thus see that increasing the Coulomb interaction U 
the HOMO and LUMO levels approach each other, in 
agreement with the general result of Section |IIB| 

The density at the terminal site is unity if the molecule 
is charge neutral. However, since the molecule with one 
electron more/less induces an impurity-like potential ^U, 
the terminal site density in this case changes according 
to 
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Appendix A: Other exact results of the Image 
Charge Model in the uncontacted case 



ni{u) 



Aii{uj,u)duj, 



(Al) 



where Aii{uj,u) ~ — ^lm[Gfi{uj , u)] is the spectral func- 
tion projected on the terminal site with an impurity-like 
potential u = ^U. The Green's function can be calcu- 
lated explicitly from the Dyson equation and reads 



Gfi(a.,C/) = G?iV)/(l-C/G;i^M) 



(A2) 



In this Appendix we derive some simple analytic result 
for the ICM with A' = A'^ = 0, C/' = and C/'^ = U. We 
will show that the main qualitative features of the sys- 
tem in equilibrium can be captured already by consider- 
ing leads of finite length. Let us consider the molecule 
with an extra electron on the LUMO level and a right 
lead with only two sites. The extra electron induces an 
impurity-like potential U on the terminal site and the 
single-particle eigenvalues of the lead Hamiltonian are 
then given by £1^2 = T\/ (U /2)2 + b'^. Let us denote by 



Here G^j^ is the unperturbed retarded Green's function 
of the semi-infinite lead and it reads 



Gil {^) 



1 

262 



iu) — sgn(w)\/w2 — 46^) ( 
(a;-iV462 -a;2) (|w| 



uj\ > 2\, 
<2\b\) 



(A3) 



If |C/| exceeds the lead hopping b a split-off state appears 
outside the energy continuum. This is illustrated in Fig. 
[25] where we plot the lead spectral function Aii{uj, U) for 
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FIG. 25: Formation of the split-off state (sharp peak in the 
bottom right panel) as U increases. In all plots b — —1.0. 



U — 0, 0.5, 1.0, 1.5. This split- ofF s tate appears as a 
pole in the Green function of Eq. (|A2[) with the energy 



e(f/) - b 



1 



, b ) 



(A4) 



Comparing the spectral structure of Fig. [25] with that 
of Fig. [20] we conclude that the extra step in the HF 
I — V curve is due to the formation of a split-ofF state 
which prevents the current to incease as the bias becomes 
larger. 



Appendix B: Density response function 

We here calculate the density response function pro- 
jected onto the terminal site of a semi-infinite TB chain 
relevant for the discussion of Section lllBI For chains with 
TVch sites the single-particle eigenfunctions and eigenen- 
ergies of the system are 4'k{i) = ( — 1)'+^ y^ ^^^^^ sin((/)fcz) 
and Efc = -26cos(0fc), where (j)k = iJ^, k = l...A'^,i 



Vch- 



By definition, the (retarded) density response function 
Xij(w) with site coordinates reads 



X.,(c.) = 2^(/,-/0 

kl 



'<P*k{i)Mi)i^k{j)Vi{3) 



(Bl) 



where, for zero temperature, fk — 0{fJ. — e^) are the 
single-particle occupations and rj is an infinitesimally 
small positive constant. Inserting the explicit form of the 
eigenfunctions and eigenvalues, changing the variables to 
k — kn /{Nch + 1) and taking the A^ch oo limit, we get 
for the i = j = 1 component 



Xii(w) 



dk dl- 
Jo ^ 



{fk - /r)_sin^fcsin^/ 
— 26(cos k — cos Z) -)- i?7 ' 



where for the half-filled system here considered /j, = 
— k). This expression can be simplified further by 
changing the variables to x = cosZ and y = cos A:. The 
integral containing /j, becomes 



dy 



-1 



V 



uj — 25(y — x) -\- ir] 



A(i)(w) -iAfiH'^), (B2) 
is the real part and 



where K^^\<^) - f dw'^-^ 



Ad) (a 



4 



dy^l-y^^\-[y-uj/{2h)f x 
e{y - uj/{2b) + 1)9(1 -{y- u/{2b))), (B3) 



is the imaginary part. Similarly one obtains the integral 

Xii containing fj. The sum xii — Xii + xii'' can now 
easily be calculated numerically. 



Appendix C: Explaining the level broadening in the 
HF approximation 

In this Appendix we show that the molecule-lead in- 
teraction in the presence of a finite hybridization renor- 
malizes the embedding self-energy already in the HF ap- 
proximation, thus explaining the broadening of the HF 
peaks in Fig. J8] Let us denote simply by G and S the re- 
tarded components of the Green function and self-energy 
respectively. For simplicity we take A' = t/' = and 
y = \, — U and we denote by 1 the terminal site 
of the right lead. We start from the Dyson equation for 
G{uj) 



{w-h- j:^^)g{oj) = 1 



(CI) 



where, in accordance with the notation of Section [II C^ h 
is the Hamiltonian in the onc-particlc Hilbcrt space and 
has the structure 



Cff hn^ 
El /iL,7 

hr,H W.L 



(C2) 



Here hr^r is the tridiagonal matrix which describes the 
right lead with matrix elements b on the upper and lower 
diagonal and zero otherwise, while hi^r, with i = L, is 
the rectangular matrix whose only non-vanishing entry 
is {hi^r)i,i = A. Projecting the Dyson equation onto HH 
and r, H we find 



(c 



{iJ-hr,r-T,f^)Gr^Hii0) 



Solving the second equation for Gr,H and inserting the 
result in the first equation we obtain the following solu- 
tion for Ghh 



GHHi(^) 



1 



-'HH 



(C3) 
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For the equilibrium system we can always choose the 
HF orbital to be real valued and therefore Im[I]^^] = 0. 
Then, inserting Eq. (|C6|) into Eq. (|C5|) and solving for 



S"! we find 



where 



VHF 



UC 



(C7) 



FIG. 26: The value of the renormahzation constant as a 
function of U for b — —1.0 and A = —0.2. The inset shows 
the dependence of the factor C as a function of U. 



C = 2 



duj _^ 
2^ 



-Im 



G'ii(a;)G_H-if(w) 



(C8) 



with 

J:„h - ^Wh + (A + I]g^)Gii(^)(A + Sfl). (C4) 

In the above equation Gn is the (1,1) matrix element of 
the Green's function of the uncontacted system (A = 0) 
with the same HF self-energy, i.e., Gr,r — ^/{'^ ~ hp ,. — 
S^^). Note that the only non- vanishing entry of the self- 
energy in the lead is — ■ Next we observe 
that the nonlocal HF self-energy can be written as 



did 

2n 



(-2zIm[GiHM]), 



(C5) 



and similarly for 'S^h- Fi^o™ the projected Dyson equa- 
tion we have 



Gm(w) = GiiH(A + Sf^)GHH(c^ 



(C6) 



This result together with its analogous for S]^^[ allows us 
to cast the self-energy in Eq. ( C4 ) in the form 



^HF 



UC 



(C9) 



where T,'^^j^{u!) = A^Gii(w) is the embedding self-energy 
of the non-interacting system. Thus, the molecule-lead 
interaction renormalizes the embedding self-energy and 
increases the broadening of the HF spectral peaks. The 
value of the constant G in Eq.(C8l can be determined 
numerically. In Fig. 26 we display (1 — CU)^"^ and G 
(inset) as a function of U . We see that G is roughly 
constant for small U . 
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